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Abstract. Let V be a periodic potential on IE 3 that is smooth everywhere except at a discrete 
set S of points, where it has singularities of the form Z/p 2 , with p(x) = \x — p\ for x close 
to p and Z is continuous, Z(p) > —1/4 for p £ S. We also assume that p and Z are smooth 
outside cS and Z is smooth in polar coordinates around each singular point. Let us denote 
by A the periodicity lattice and set T := R 3 /A. In the first paper of this series 1201 . we 
obtained regularity results in weighted Sobolev space for the eigenfunctions of the Schrodinger- 
type operator H = — A + V acting on L 2 (T), as well as for the induced k-Hamiltonians 
obtained by resticting the action of H to Bloch waves. In this paper we present two related 
applications: one to the Finite Element approximation of the solution of (L + H^)v = f 
and one to the numerical approximation of the eigenvalues, A, and eigenfunctions, u, of H^. 
We give optimal, higher order convergence results for approximation spaces defined piecewise 
polynomials. Our numerical tests are in good agreement with the theoretical results. 
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1. Introduction and statement of main results 

In this paper, which is the second part of work begun in |20) . we present applications of the the- 
oretical regularity results in the first part of this paper to Finite Element Method approximation 
schemes. The first application is to approximation of eigenvalues, A, and eigenfunctions, u, of the 
Bloch operator, H^, associated to a periodic Hamiltonian operator with inverse square potential 
at isolated points. For example, one of our main results, Theorem |1.1| yields optimal orders of 
convergence for the Finite Element approximations of the eigenvalues of using graded meshes. 
These rates are higher than those that can be obtained using standard meshes. The second ap- 
plication is to the Finite Element Method, again using graded meshes, applied to equations of 
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the form (L + H-^)v = /. The final section of this paper presents numerical tests showing good 
agreement with our theoretical results for this second problem. 

Hamiltonian operators with inverse square potentials arise in a variety of interesting contexts. 
The standard example of a Schrodinger operator with c/p potential is a special case of the inverse 
square potentials wc consider, where the function p 2 V vanishes to order 1 at the singularity, and 
the results of this work apply to such operators. But in addition, Hamiltonians with true inverse 
square potentials arise in relativistic quantum mechanics from the square of the Dirac operator 
coupled with an interaction potential, and they arise in the interaction of a polar molecule with 
an electron. See [29j |32] for further applications of inverse square potentials to physics. See also 
[UUHl [25, 28, 9] for related results on operators with singular coefficients. Thus it is interesting in 
several areas of physics to understand how to approximate solutions to equations involving such 
operators. 

Before we can state our approximation results, we must fix some notation and state the assump- 
tions we make about our Hamiltonian operators. Consider a Hamiltonian operator H := —A + V 
that is periodic on R 3 with triclinic periodicity lattice A. Its fundamental domain is a parallelop- 
iped whose faces can be identified under the symmetries of H to form the torus T = K 3 /A, which 
is how we will denote this fundamental domain in the remainder of this paper. Let p(x) be a 
continuous function on T that is given by p(x) — \x — p\ for x close to p, is smooth except at the 
points of S, and may be assumed to be equal to one outside a neighbourhood of S. 

We need two assumptions about the potentials V that we will consider in this paper. First, we 
assume that V is smooth except at a set of points S C T, near which it has singularities of the 
form Z/p 2 , where Z is continuous on T and smooth in polar coordinates around p. We denote 
this as follows. 

(1) Assumption 1 : Z := p 2 V G C(T) n C°°(T^5). 

Assumption 1, more precisely the continuity of Z at S, allows us to formulate our second assump- 
tion. Namely, 

(2) Assumption 2 : rj := min yI/4 + Z(p) > 0. 

pes 

In particular, we assume that for all p G S, Z(p) > —1/4. These assumptions are sharp in the 
sense that the analysis yields fundamentally different results if either one fails. In particular, the 
value rj = —1/4 corresponds to the critical coupling for an isolated inverse square potential in 
K 3 where the system undergoes a transition between the conformal and non-conformal regimes 
[2"5] , If the first assumption fails, then the available analytic techniques are much weaker, see for 
instance |171I16] . In either case, the approximation theorems in this paper fail if either assumption 
is violated. More details of this are included in the first part, [221, and a study of the analysis 
when these assumptions are relaxed will be examined in a forthcoming paper. 

We are interested in understanding the spectrum and generalised eigenfunctions of the operator 
H. As usual, we do this by studying Bloch waves. Recall that if k is an element of the first Brillion 
zone of A, that is, is an element of the fundamental domain of the dual lattice of A, then a Bloch 
wave with wave vector k is a function in L ; 2 oc (IR 3 ) that satisfies the semi-periodicity condition 

(3) Vk(z + X) = e ik ' x -0 k (» VA G A. 
It is well known that such a Bloch wave can be written as 

(4) Vk(z) = e^ x u k (x) 

for a function that is truly periodic with respect to A and thus can be considered as living on 
the three-torus T. We define the k-Hamiltonian on L 2 (T) by 



(5) 



3 

H k :=-J2(d j +ik j ) 2 + V. 

3=1 
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Then we have further that if a Bloch wave V'k is a generalized eigenfunction of H with generalized 
eigenvalue A, then the function itk := e^ l]i ' x ip^(x) is a standard L 2 -eigenfunction of with 
eigenvalue A. Let Xj, j > 1, be the eigenvalues of ifk, arranged in increasing order, .. . < Xj < 
Xj+i < . . ., and repeated according to their multiplicities. That is, if E(X) denotes the eigenspace 
of H k corresponding to A, then A is repeated dim(S(A)) times. 

As usual, for our finite clement approximation results, we consider a sequence S n of finite 
dimensional subspaces of the domain of i?k and let R n denote the Riesz projection onto S n , that 
is, the projection in the bilinear form ((L + H-^y, w)l 2 (t-^s)> (for a suitable C > 0), associated 
to H^. Let -f/k.n := RnHkR n be the associated finite element approximation of H^, acting on 
S n . Denote by Xj t „ the eigenvalues of the approximation H^ n , again arranged in increasing 
order, . . . < Aj.„ < Aj +1 „ < . . ., and repeated according to their multiplicities and let Uj iU G S n 
be a choice of corresponding eigenfunctions (linearly independent). The spaces S n we use for 
our theorems are defined in terms of a sequence of graded tetrahedral meshes T n := k n (7~o) 
on T (sometimes called triangulations) , given by sequential refinements, associated to a scaling 
parameter k, of an original tetrahedral mesh 7o- We describe the meshing refinement procedure 
in detail in Section |3j We will take S n = S(T n ,m), the finite element spaces associated to these 
meshes (i.e., using continuous, piecewise polynomials of degree m). 

Our first theorem, which is a theoretical result for the finite clement method approximation of 
eigenvalues and eigenfunctions of using tetrahedralisations with graded meshes, is as follows. 

Theorem 1.1. Let Xj be an eigenvalue of and fix < a < f], a < m. Let Aj,„ be the finite 
element approximations of Xj associated to the nested sequence T n of meshes on T defined by 
the scaling parameter k = 2~" l / a and piecewise polynomials of degree m. Also, let Uj_ n be an 
eigenbasis corresponding to Xj^ n . Then there exists a constant c(Xj , a) independent of n such that 
the following inequalities hold for a suitable eigenvector Uj € E(Xj): 

|Aj-A J>B | <c(A J ,a)dim(5„)- 2m/3 , 

11% - M j,nllK;}(T^5) < c ( A j> a ) dim(S'„) _m/3 , 
where the space K,\(T \ S) is a weighted Sobolev space defined below in Equation^ 

For our second theorem, we consider the finite element approximations of the equation 

(6) (L + H k )v = f, forL>C , 



where Co is the constant from Theorem 2.1 below. We then define the form a(y,w) :— ({L + 
Hk)y,w) and let v be the solution of Equation ([6| above. We then define the usual Galerkin 
Finite Element approximation v n of v as the unique «„ £ S„ := S(T n , m) such that 

(7) a(v n ,w n ) := ((L + H k )v n ,w n ) = (f,w), for all w n <E S n . 

Then Theorem |l.l| ;ogether with the Lax-Milgram Lemma and Cea's lemma imply that we 
have the following h m quasi-optimal rate of convergence. 

Theorem 1.2. The sequence T n '■= k n (To) of meshes on P defined using the k-refinement, for 
k = 2~ m / a , < a < r\, a < m, and piecewise polynomials of degree m, has the following 
property. The sequence v n € S n := S(%i,rn) of Finite Element (Galerkin) approximations of v 
from Equation ^ satisfies 

(8) \\v - v„||k;;oks) < Cdim(5„)- m / 3 ||/||^ ril(T ^ ) , 

where C is independent of n and f. 

These theorems are interesting because it is known that the convergence rate of a standard 
finite element method (i.e., based on quasi-uniform meshes) is limited. However, under the 
assumptions on our potentials, if we use graded meshes instead, we can obtain an approximation 
rate as fast as we like by using polynomials of sufficiently high degree in the elements. This is due 
to the fact that although regularity of the associated Bloch waves is limited in terms of standard 
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Sobolev spaces on T, it is arbitrarily good with respect to weighted Sobolev spaces. We will recall 
the definition of these spaces and the relevant regularity results from [50] along with additional 
background in Section [2] 

The remainder of the paper is organised as follows. In Section [3] we first describe the k- 
refinement algorithm for the three dimensional tctrahcdral meshes, which results in a sequence 
of meshes T n - We then prove a general interpolation approximation result for the sequence of 
finite element spaces associated to this sequence of meshes. In Section [4] we use this general 
approximation result to prove our main approximation results This section includes in particular 
the proofs of Theorem |1.1| and Theorem L2 as well as an additional result about the condition 



number of the stiffness matrix associated to the finite clement spaces, S n . In the last section, 
Section [5j we discuss results of numerical tests of the method for solving equations of the form 
(L + Hk)v = f and compare them to the theoretical results. 

Acknowledgements. We would like to thank Bernd Ammann, Douglas Arnold, and Catarina 
Carvalho for useful discussions. We also thank the Leverhulme Trust whose funding supported 
the fourth author during this project. This project was started while Hunsicker and Nistor were 
visiting the Max Planck Institute for Mathematics in Bonn, Germany, and we are greatful for its 
support. 



2. Background results 

In this section we recall some definitions and results from 20 , as well as the classical approx- 
imation result for Lagrange interpolants (see [TTJ [TSl 131])) that will be used in the proofs of 
the approximation theorems above. First, these results arc given in terms of weighted Sobolev 
spaces which are defined as follows: 

(9) /C™(T n S) := {v : T \ S -> C, pW-'^v G £ 2 (T), V \/3\ < to}. 

These spaces have been considered in many other papers, most notably in Kondratiev's ground- 
breaking paper [24] . 

The first result that we recall guarantees the existence of solutions of equations of the form 
[L + H\t)v = f for L greater than some constant Co, and identifies the natural domain of H^. Let 
us fix smooth functions \ p supported near points of S such that the functions \p have disjoint 
supports and x p = 1 in a small neighbourhood of p G S. Then Theorem 1.1 and Proposition 3.6 
from |20j combine to give right away the following result. 

Theorem 2.1. Let V be a potential satisfying both Assumptions 1 and 2. Then there exists Co > 
such that L + H\s_ : /C^A (T \ S) — > KJ^Z-y (T % S) is an isomorphism for all m G Z>o, all \a\ < n, 
and all L> C* . Moreover, for any u G JC^(T \ S) satisfying (L + H k )v = f G H m - l {T \ S), 
we can find constants a p G R such that 

u reg :=u-J2 XvP^ 11 ^- 1 ' 2 G v S). 

pes 

We obtain, in particular, that has a natural self-adjoint extension, the Friedrichs extension. 
Therefore, from now on, we shall extend to the domain of the Friedrichs extension of L + H^, 



as in the above Theorem. Let us denote by V{H k ) its domain. Then Theorem 2.1 gives that 
f(-ffk) = /C| (T \ S) for min p Z(p) > 3/4, and, in general, 

(10) P(ff k )c/CL(T\5), for a < r) := min Jill + Z{p) and a < 1 

p 

so that I?(iJk) C K,\(T \ S) C i7 1 (T \ S), since we assumed that min p Z(p) > —1/4. 

We can now state a regularity theorem for the eigenfunctions of near a point p G S, or 
equivalently, for Bloch waves associated to the wavevector k. 
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Theorem 2.2. Assume that V satisfies Assumptions 1 and 2 and let u G 2?(£/k) satisfy H^u = 
Xu, for some Aef. Then we can find constants a p G R such that 

u-Y x P P^ 1/i+Z(p) - 1/2 G K™X\{T \ S), Va' < minV9/4 + %). 

z — ' p£5 

in particular, u G /C^\ 1 (T \ 5), where a < r\ :— min pe 5 yl/4-j- Z(p) and m G Z + is arbitrary. 



See also [22, 23 for some related classical results in this area. Theorems 2.1 and 2.2 lead to 
an estimate for the distance from an element in the domain of i?k to the approximation spaces 
that we construct using graded meshes. 

Next, recall the definition of Lagrange interpolants associated to a mesh. Let us choose P to 
be a parallelopiped that is a fundamental domain of the Lattice A. That is, K 3 = U ?/6 a(2/ + P) 
and all y + P disjoint. Let T = {T^} be a mesh on P, that is a mesh of P with tetrahedra T,. 
We can identify this T with a mesh T' of the fundamental region of the lattice C (that is, to the 
Brillouin zone of C). Fix an integer m G N that will play the role of the order of approximation. 
We denote by 5(7", m) the finite element space associated to the degree m Lagrange tetrahedron. 
That is, 5(7", m) consists of all continuous functions x : P — >■ K such that x coincides with a 
polynomial of degree < m on each tetrahedron TeT and \ is periodic. This means the values 
of x °n corresponding faces coincide, so x will have a continuous, periodic extension to the whole 
space, or alternatively, can be thought of as a continuous function on T ■ We shall denote by 
wj = wix G 5(7", m) the Lagrange intcrpolant of w 6 H 2 (M. 3 ). Let us recall the definition of 
wjj-. First, given a tetrahedron T, let [to, t\, ti, t 3 ] be the barycentric coordinates on T. The 
nodes of the degree m Lagrange tetrahedron T are the points of T whose barycentric coordinates 
[to, ti, t2, tg] satisfy mtj € Z. The degree m Lagrange interpolant Wi,t of u is the unique function 
Wip- £ 5(T, m) such that w = wij- at the nodes of each tetrahedron T G 7~. The shorter 
notation Wj will be used when only one mesh is understood in the discussion. 

The classical approximation result for Lagrange interpolants ( [H [H] [T51 [21] ) can now be stated. 

Theorem 2.3. Let T be a mesh of a polyhedral domain P C R 3 with the property that all 
tetrahedra comprising 7" have angles > a and edges < h. Then there exists a constant C{ct, m) > 
such that, for any u G iJ m+1 (P), 

II" - u/IIhHp) - C{a,m)h m \\u\\ Hm+ i {r) . 

Finally, we recall two properties of functions in the weighted Sobolev spaces /C™ (T \ S) that are 
useful for the analysis of the approximation scheme we use with graded meshes. The proofs of these 
lemmas are contained in |21) and are based on the definitions and straightforward calculations. 

Lemma 2.4. Let D be a small neighborhood of a point p G S such that on D, p is given by 
distance to p. Let < 7 < 1 and denote by "fD the region obtained by radially shrinking around 
p by a factor of 7. Then 

IMIiC2»(i?) = (j) a ~ 3/2 \H\K?(iD)- 
Lemma 2.5. If m> m' , a> a' and < p < 8 on D, then 

\W\ K7 ; {D) <6 a - a '\\w\\^ {D) . 

We can now continue to the definition of the mesh refinement technique and the proof of the 
general approximation theorem underlying our two main theorems. 

3. Approximation and mesh refinement 

Our two main theorems follow from standard results, such as Cea's Lemma (for the proof of 
Theorem 1.2 1 and the results used in [31 [51 |H [101 [3D] (for the proof of Theorem |1.1[ ), together 



with the following underlying approximation theorem: 
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Figure 1. The initial tetrahedron {xq, x\, x%, £3} (left); eight sub-tetrahedra 
after one fc-refinement (right), k = ^ = « = ^f. 

Theorem 3.1. There exists a sequence T n of meshes ofT that depends only on the choice of a 
parameter k < 2~ m / a , < a < r\ and a < m, with the following property. If u G /C^ L 1 (T \ S), 
then the modified Lagrange interpolant uij- n € S(7~ n ,m) of u satisfies 

\\ u - u LtJk\(T^S) < C7dim (' S 'n)~ m/3 |l'"llK™ + + 1 (TxS)' 

where C depends only on m and a (so it is independent of n and u). 

In this section we will define the mesh refinement process and prove Theorem |3.1| The first step 
is to describe the refinement procedure that results in our sequence of meshes (or triangulations) . 
This is based on the construction in [B] and in [5], thus we refer the reader to those papers for 
details, and here give only an outline and state the critical properties. The second step is to 
prove a sequence of simple lemmas used in the estimates. The third step is to prove the estimate 
separately on smaller regions. This uses the scaling properties of the meshes in Lemmas |2.4| and 
[23] together with Theorem [231 

3.1. Construction of the meshes. We continue to keep the approximation degree m fixed 
throughout this section. Fix a parameter a and let k = 2~ m / a . In our estimates, we will chose 
a such that a < r\ :— min p ^/l/4 + Z[p) and a < m. Let I denote the smallest distance between 
the points in S. Choose an initial mesh 7o of P with tetrahedra such that all singular points of V 
(i.e., all points of S) are among the vertices of 7o and no tetrahedron has more than one vertex 
in S. We assume that this mesh is such that if Fi and F 2 are two opposite faces of P, which 
hence correspond to each other through periodicity, then the resulting triangulations of F\ and 
i*2 will also correspond to each other, that is, they are congruent in an obvious sense. 

We start with a special refinement of an arbitrary tetrahedron T that has one of the vertices 
in the set S. Our assumptions then guarantee that all the other vertices of T will not be in S. 
Motivated by the refinement in [51151 112L 1131 127j . we define our k-refinement algorithm for a single 
tetrahedron that divides T into eight sub-tetrahedra as follows. 

Algorithm 3.2. ^-refinement for a single tetrahedron: Let Xq, x\, x%, x% be the four vertices 
of T. 

We denote T by its vertex set {xo,xi,X2,Xa}. Suppose that xq g 5, so that xq is the one 
and only vertex that will be refined with a ratio k e (0,1/2]. We first generate new nodes 
Xij, < i < j < 3, on each edge of T, such that Xij = (Xi + Xj)/2 for 1 < i < j < 3 and 
x j = (1 — k)x + kxj for 1 < j < 3. Note that the node Xij is on the edge connecting X{ and Xj. 
Connecting these nodes Xij on all the faces, we obtain 4 sub-tetrahedra and one octahedron. The 
octahedron then is cut into four tetrahedra using x\% as the common vertex. Therefore, after one 
refinement, we obtain eight sub-tetrahedra (Figure [T]), namely, 

{£0,2:01, £02, £03}, {xi, £01, £12, £13}: {x2,x 2,x 12 ,X2i\, {^3,a;o3,a;i3,a;23} 

{xqi, X 2, X 03 , X13}, \Xqi, X 2, X\2i £13}, {^02, ^03, ^13, ^23}, {^02, ^12, ^13, X23}. 
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Algorithm 3.3. fc-refinement for a mesh: Let T be a triangulation of the domain P such that 
all points in S are among the vertices of T and no tetrahedron contains more than one point 
in S among its vertices. Then we divide each tetrahedron T of T that has a vertex in S using 
the fc-refinement and we divide each tetrahedron in T that has no vertices in S using the 1/2- 
refinement. The resulting mesh will be denoted fc(T). We then define T n = k n (%), where 7o is 
the intial mesh of P. 

Remark 3.4. According to [8 j, when fc = 1/2, which is the case when the tetrahedron under 
consideration is away from iS, the recursive application of Algorithm |3.2| on the tetrahedron 
generates tetrahedra within at most three similarity classes. On the other hand, if fc < 1/2, 
the eight sub-tetrahedra of T are not necessarily similar. Thus, with one fc-refinement, the 
sub-tetrahedra of T may belong to at most eight similarity classes. Note that the first sub- 
tetrahedra in Algorithm |3.2| is similar to the original tetrahedron T with the vertex xo £ S and 
therefore, a further fc-refinement on this sub-tetrahedron will generate eight children tetrahedra 
within the same eight similarity classes as sub-tetrahedra of T. Hence, successive fc-refinements 
of a tetrahedron T in the initial triangulation 7o will generate tetrahedra within at most three 
similarity classes if T has no vertex in S. On the other hand, successive fc-refinements of a 
tetrahedron T in the initial triangulation will generate tetrahedra within at most 1 + 7 x 3 = 22 
similarity classes if T has a point in S as a vertex. Thus, our fc-refinement is conforming and 
yields only non-degenerate tetrahedra, all of which will belong to only finitely many similarity 
classes. 

Remark 3.5. Recall that our initial mesh 7o has matching restrictions to corresponding faces. 
Since the singular points in S are not on the boundary of P, the refinement on opposite boundary 
faces of P is obtained by the usual mid-point decomposition. Therefore, the same matching 
property will be inherited by T n - In particular, we can extend T n to a mesh in the whole space 
by periodicity. We will, however, not make use of this periodic mesh on the whole space. 

For each point p G S and each j, we denote by V P j the union of all tetrahedra of 7} that have 
p as a vertex. Thus V P j is obtained by scaling the tetrahedra in V p o by a factor of k J with center 
p. In particular, the level n > j refinements of 7o give rise to a mesh on lZ p j := V p (j_i) \ V P j. 
Define 

0:=P\U pe5 V p o- 



According to Definition 3.3 both and Ll p ^sV p o are triangulated using the fc-refinement. For 
each tetrahedron in V p o we use the fc-refinement for a single tetrahedron, while for fi we use the 
1/2-refinement for meshes, which is, of course, a uniform refinement. Then, we can decompose 
P as the union 



(11) P = U pes ^ U? =1 U pj U V pn 

where each set in the union is a union of tetrahedra in T n , 

Remark 3.6. Note that the size of each simplex of T n contained in 51 is C(2 _n ), the size of each 
simplex of T n contained in TZ P j is 0(k^2~^ n "^), and the size of V pn is 0(k n ). In addition, the 



number of tetrahedra in T n is C(2 3 ™) (see Algorithm 3.3) 



We now define the finite element approximation u n G S(T n , m) to the equation (L + H^jv = f, 
where 7~ n is obtained by applying n times the fc-refinements to 7o, where fc = 2 _m / a , < a < rj, 



a < m, and A > satisfies Theorem 2.1 Then u n is defined for any v n € S(T n ,m) by 



(12) (H k u n ,v n ) +X(u n ,v n ) := (Vit„, Vv n ) L 2 + ((V + X)u n , v n ) L i = (f,v n ) L 2. 



Note that Theorem 2.1 gives that the finite element solution u n € S(T n , m) C JC\ is well defined 
by (12). The approximation properties of u n are discussed in Theorem 



8 



E. HUNSICKER, H. LI, V. NISTOR, AND V. USKI 



3.2. Proof of Theorem |3.1| , Note that the singular expansion of Theorem |2.2| shows that the 
value of an eigenvalue u of at a singular point in S may not be defined. Therefore, we must 
define the modified degree m Lagrange "interpolant" ui )Tl = u/,t„ associated to the mesh T n , 
such that 

Ui, n {x) = u(x) for any node x S 
Ui lV ,(x) = if X G S. 



(13) 



Alternatively, we can take the modified Lagrange interpolant to be zero on the whole tetrahedron 
that contains a singular point. By construction, the restriction of T n to lZ P j scales to the restriction 
of Tn-j+i to 1Z P \. From now on, we refer to uj n — uj j- n as the modified interpolation defined 
in (13 1. The following lemma is based on the definition of the fc-refinement and the discussion in 
Remark 13.61 

Lemma 3.7. For all x G TZpj, Ui jn {x) = u/ )ri _j+i(fc - W _1 '(i)), where k~^~^(x) := p + (x — 
p)/kSi~ l > is the dilation with ratio k~^~^ and center p. 



Recall that p 2 V G C°°(T \ S) nC(T) and min,, Z{p) > -1/4. That is, V satisfies Assumptions 
1 and 2. 

We can now give the proof of Theorem |3.1[ 

Proof. Recall that V p o consists of the tetrahedra of the initial mesh 7o that have p as a vertex 
and that all the regions V p are away from each other (they are closed and disjoint). We used this 
to define fi :=P\ U p V p o- The region V p j is obtained by dilating V p with the ratio k? < 1 and 
center p. Finally, recall that lZ P j — V p (j-\) \ V P j. Let R be any of the regions fl, lZ P j, or V pn . 
Since the union of these regions is P, it is enough to prove that 

h-ui,T n \\jcl(R^s) < C d™(Sn)~ m/3 \\u\\ K rn.+i( R ^ s) , 

for a constant C independent of R and n. The result will follow by squaring all these inequalities 
and adding them up. In fact, since dim(5„) _m/ ' 3 = 0(2~™ m ), it is enough to prove 



(14) \W-u IiT Jki(R-,s)<C2- 
again for a constant C independent of R and n. 



if r = n 



\ UpVpO; the estimate in (14) follows right away from Theorem 2.3 



For the 



other estimates, recall that < k < 2 m / a : where < a < r\ and a < m. We next establish the 
desired interpolation estimate on the region R — 7Z P j, for any fixed p G S and j = 1, 2, . . . , n. 



Let u(x) — u(k 3 x). From Lemmas 2.4 and 



3.7 



h-^,»ki( RM .) = (^Y /2 ||£- 



Since /C" l (7\t p i) is equivalent to H m {1Z p i), we can apply Theorem 
get 

(15) Ik-ti/.nkK**) < C{k i ~ l ) l/2 2~ m{ - n ~ i+r> \\u\\ 

Now applying Lemma 
hand side in ( 15 ) 



we have 

(ui,n)\\lC{(K pl ) = ( fcJ ~V /2 |l«- W/,„_ j + i||k;1(TC p1 )- 

with h = 0(2-("-J' +1 )) to 



2.3 



2.4 



to scale back again and using also k — 2 m / a , we get that the right 
C(F _1 ) 1/2 2 _m( "^ :,+1) ||u|| K;m +i(K pl ) = C{k j - x ) a 2- m ^ n - j+l \ 



^a+1 C^pj) 

< C2- 



I K) ■ 



This proves the estimate in ( 14 ) for R = 1Z 

It remains to prove this estimate for R — V pn . For any function w on V pn , we let w(x) — w(k n x) 
be a function on V p . Therefore, by Lemma 2.4 



(16) 



\\ u - u I,n\\K{(V pn ) = { kn ) 1/2 \\ u ' rr uiA\K.\(V p ) 
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and by 3.7 (which follows from the definition of the meshes 7~k and from the fact that interpolation 

(/c n ) 1/2 l 



commutes with changes of variables) 
(17) (fc™) 1/2 || 



Ul,n\\lCl(V p ) 



ui,o\\kI(v p )- 



Now let x be a smooth cutoff function on V p such that x = in a neighborhood of p and = 1 at 
every other node of V p . 



Define v := u — x u - Then, by (13), 
(k n f/ 2 \ 



U IMK{(V P ) 



< 



(18) 



(fc n ) 1/a ||* + Xfl-fi/,o||x:i(v,) 
(k n ) 1/2 (\\v\\ KliVp) + \\xu- ui,o\\>cl(y p )) 
(k n ) 1/2 (\\v\\ KliVp) + \\xu- (xv)i,o\\k\(v p ))- 



Since x vanishes in the neighborhood of p we can consider multiplication by x as P°° times a 
degree b-operator. Thus it is a bounded operator on any weighted Sobolev space. Thus 

(19) \\v\\k\(v p ) < INkf (V„) < \\u\\k.?(v p ) + llx^lkf (V p ) < C\\u\\ KT{Vp) , 

where C depends on m and, through %, the nodes in the triangulation. 



Using (|16j), ([17j), (|18j), (|19j), Lemma |2.5[ and Theorem [^3) we have 

C(fc™) 1/2 



u I,n\\K\(V p7l ) 



< 
< 
< 
< 



^k;(v p ) + - {Xu)i,o\\ K \{v p ) 
C{k n ) 1/2 {\\u\\K\(v P ) + \\xu\\ H ^(v p) ) 



\V P )> 



C(fc") i/2 (||«k; ( v p ) + NI, 



< C2- 



This proves the estimate of Equation ( 14 ) for R = V pn and completes the proof of Theorem 

□ 



4. Applications to Finite Element Methods 

We can now turn to the proofs of the theorems stated in the introduction. First, Theorem 
1 1.1 1 follows from our general approximation result, Theorem 3.1 and the standard results on 
approximations of eigenvalues and eigenvectors (eigenfunctions in our case) discussed, for instance, 
in [5J [SI HI [TUl [30] . More precisely, using the notation introduced in the introduction, we have the 
following. Let us denote by E(X) the eigenspace of corresponding to the eigenvalue A and by 
-Ei (A) C E(X), the subspace consisting of vectors of length one. Then the following result is well 
known (see for instance Equations (1.1) and (1.2) in [5]). We state it only for our operator H^, 
although it is valid for more general self-adjoint operators with compact resolvent. 

Theorem 4.1. There exists a constant C > with the following property. Let V C /Cj(T \ S) 
be a finite dimensional subspace and R : /Cj (TT \<S) — > V the projection in the energy norm. Let 
Wj t n € V be an eigenbasis of RH^R, namely RH-^Rwj^n ~ RH-^Wj^ ~ Xj, n uij, n , with the \j ;n 
arranged in increasing order in j . Then 



and, furthermore 



\Xj - \j, n \ < C sup inf \\u - xIUoks) 



fj.nll/cirTvS) < c SU P inf Hu-xll/cKT-vS), 



for a suitable eigenvector Uj £ E(Xj). 

The proof of Theorem |1.1| will then be obtained from Theorem |4.1| as follows. 
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Proof, (of Theorem 1.1 1. We need to estimate sup ueEi ( X ) ^ n ^xeS„ \\ u ~ x\\jc{- To this end, let 

(/j + X)u. Theorem 



2.1 



then 



us notice that any u £ E(X) C K.\(T \ S) satisfies (// + H^)u 

gives llull^m+i < Cm^llwll^m-i for a suitably large fi that depends on A and a < rj. A bootstrap 
argument then gives for any u £ E(X) that HuH^-r-, 
u € Ei{\j) (thus ||w||/c; = 1), the following. 



n < C: 



m.A I 



\K\- 



Theorem 



3.1 



then gives for 



sup inf ||«-xIIki(Tv.s) < SU P \\ u ~ u i.tA\k\(i^s) 



< C sup dim(S f „) 
The proof of Theorem |1.1| is now complete. 



-m/3| 



I /C^+j 1 (T 



„s) < c(m, Aj) dim(S„) 



-m/3 



□ 



Next, the proof of Theorem 1.2 follows from Theorem 3.1 the Lax-Milgram Lemma and Cea's 
lemma. We note some consequences of this theorem. 

Remark 4.2. First, in the case / € i/ m_1 (T \ <S>), by the estimate in Equation p]), we have 

Il«-««llx:icrv5) < Cdim^n)—/ 3 !!/!!^^^ < Cdim(5 n )- m / 3 ||/|U m - 1(TxS)) 
as long as the index in Theorem |1.2| is chosen such that < a < 1 . 



As in the classical Finite Element Method, a duality argument yields the following L 2 -convergence 
result. 



Theorem 4.3. In addition to the assumptions and notation in Theorem assume that < 
a < 1. Then the following L 2 estimate holds 

\\v-v n \\ LHT) < Cdim(S„) ( - m - 1)/3 ||/|| ffn 



Proof. We sketch the proof by using the duality argument in weighted Sobolev spaces. Consider 
the equation 



(20) 



(L + H k )w 



in T. 



(So we use periodic boundary conditions on P.) The definition of the Galerkin projection v n of 
v, Equation Q, then gives 

(v -v n ,v- v n ) = ((L + Hic)w, v - v n ) = {{L + H]t)(w - w n ), v - v n ), 

is the finite elemen 

C\\v - w„||z, 2 (T) by Theorem 



where w n is the finite element solution of Equation ( |20[ ) on T n , We also have ||w||yc 2 t (TxS) ^ 

v n e L 2 {T) C /C°_ X (T \ S). Therefore, applying 



Theorem 



1.2 



2.1 



since !' 



to v 



€ L 2 (T) and m = 1, we have 



\\v - V n \\ L 2 {Tj < C\\w-w n \\,ci m \\v-v n \\ K i m /\\v-v n \\ L 2(T) 



< Cdim(S* r 



\-V3| 



^nlkj(T) < C dim(5' n ) 



(-m-l)/3i 



-i(T)- 



This completes the proof. 



□ 



4.1. Condition number of the stiffness matrix. It is important that the discrete system S n 
that we use is well-conditioned for us to be able to realise the theoretical approximation bounds 
in practice. Thus we need to additionally obtain upper and lower bounds on the eigenvalues of 
the stiffness matrix that arises in calculation. 

Recall the standard nodal basis function <j)j of the space S n := S(T n , m). It consists of functions 
that are equal to 1 at one node and equal to zero at all the other nodes. For convenience, we now 

— 1/2 

instead consider the rescaled bases ipj := hj <pj, where hj is the diameter of the support patch 
for <j)j. Then, we consider the scaled stiffness matrix 

(21) 



An := (a(ipi,ipj)) 
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from our graded finite element discretization ([T]). In practice, A n can be obtained from the usual 
stiffness matrix (a(<f>i, (f)jj) by a diagonal preconditioning process. We point out that similar scaled 
matrices were considered in I26| for condition numbers of other Galerkin-based methods. 

For a symmetric matrix A, we shall denote by X m ax{A) the largest eigenvalue of A and by 
\ min (A) the smallest eigenvalue of A. Thus the spectrum of A is contained in [X min (A), X max (A)], 
but is not contained in any smaller interval. We first have the following estimates regarding 
properties of functions. 

Lemma 4.4. Let Tj be a tetrahedron in the mesh T n and let diam(Tj) denote the diameter ofTi. 
Then, for any (x n S S n and fi g there exists a constant C > independent of n, \i n and 

fx, such that 

(22) IKIIffi(T 4 ) < CdiamtTO^HMnlU-^) < C\\ 

(23) IImIU 6 (o) < C||^||iji(n). 

Furthermore, writing /x n = ^Cjtpj, where (fj :— h- 1 ^ 2 (j}j are rescaled basis functions, then 

(24) C- 1 ' 2 c ) ^ diam^OIKIli-m) < c E c l 

je.node{Ti) j£node(Ti) 



Proof. We shall show ( 22 1 and ( 24 1 since (|23[) is a standard result in [TS] . Recall that all the 



tetrahedra Tj belong to a finite class of shapes (or similarity classes) in our graded triangulation. 



Thus, the bounded constant C in (22 1 follows from the inverse estimates in [TT] I14j. 



As for (24), note fj, n = ^Ci^pi — J^Ci^i- Based on the definition of the basis function ipt and 



of the graded mesh, 



(25) C- 1 diam(T l ) 4 1/2 c l < a < Cdiam(T l )- /2 

On the reference tetrahedron T, both and (X)jgnode(T) ^j) 1 ^ 2 ar e norms for the finite 

element function v\f, where v is obtained by the usual scaling process and the summation on Cj 
is for all the nodes in T. Based on equivalence of all norms for a finite dimensional space, we have 

c( £ < \\v\\ Loe{t) < c{ J2 ^) 1/2 - 

jSraode(T) j£node(t) 



C, 



This, together with (251, implies 



C ° 2 j ^ diam ( T <)ll«lli«(r < ) < C E c l 

j£node(Ti) jenode(Ti) 

which completes the proof. □ 

Therefore, we have the following estimates on the eigenvalues of the stiffness matrix. 
Lemma 4.5. Let A n be the stiffness matrix from the finite element discretization corresponding 



to the rescaled nodal basis (pj of the space S n := S{T n ,Tn) in Equation (21). Then, 

X m ax{A n ) < M, 

where the constant M is independent of the mesh level n. 

Proof. Let us fix the mesh level n. All the constants below will be independent of n. Let {Tj} 
be the tetrahedra forming our mesh T n . Let v € S n be arbitrary and write v = • CjPj and 
V := (cj). Then, by Lemma 3.4 in [5D] we have 

V T A n V = a(v,v) < C\\v\\l m < C\\v\\ 2 H1{r) < \Mm {T& 



By the inverse inequality i22h and the estimate (24), we further have 



v T A n v < cJ2 d ™™( T *)\M 2 L°°m)<cJ2 c2 j ^ cvTv > 
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where diam(Ti) is the diameter of the tetrahedron T{. This completes the proof. □ 



Lemma 4.6. We use the same notation as the one for Lemma \4-5\ The smallest eigenvalue of 
the stiffness matrix A n , 

\ mm (A n )>Cdim(S n )- 2 / 3 . 
Proof. For any v € S n , we use the notation v = Y]j Cjipj, V := (cj), and diam(Ti) denotes the 



diameter of Tj, as before. In view of (24), the inverse estimate (22), Holder's inequality, and the 



Sobolev embedding estimate (|23|, we then have 

t> mi „.n 2 ^ \ ii.. I, 

li 6 (Ti) 



V T V = ^c^C^diam^l^ll^^^C^IHI 2 



3 * * 

i i 

< Cdim(S'„)5 \\v\\ 2 HHr) < Cdim{S n f3 V T AV. 



(P) 



□ 



Then, we have the estimate on the condition number. 



Theorem 4.7. Let A = (a(ipi,ipj)) be the stiffness matrix. Then the condition number k(A) 
satisfies 

k{A) < Cdim(S* n ) 2/3 . 
The constant C depends on the finite element space, but not on dim(S l „). 



Proof. Using k(A) = X max (A)/X m i n (A), we obtain the estimate by Lemmas 4.5 and 4.6 □ 



5. Numerical tests of the finite element method 



We now present the numerical tests for the finite element solution defined in ( 12 ) approximating 
possibly singular solutions to Equation [6j 

To be more precise, suppose that our periodicity lattice is 2Z 3 and we choose our fundamental 
domain P = [—1, l] 3 to be a cube of side length 2. We impose periodic boundary condition on 
the following model problem 

(26) (L + H k )v := (-A + ^r" 2 + L)v = 1 in Q, 

where r = |ar|, S > —1/4, L > 0, and the cut-off function tp := e r "^ r ~ r ^ +1 for r 2 < r c and 
tp = for r 2 > r c ; in the tests, we chose r c — 0.25. Note that if S > 0, it is clear that the 



operator L + is positive on K.\ (see Theorem 2.1). We use the C° linear finite element 
method on triangulations graded toward the origin with grading ratio k > (Recall that k = 0.5 
corresponds to the quasi- uniform refinement.) 

To enforce the periodic boundary condition for the finite element functions, we use meshes 
where all the boundary nodal points are symmetric about the mid-plane between opposite faces 
of the cube. Any set of the symmetric nodes will be associated to the same shape function 
in the discretization. For example, nodes on edges of the cube generally have three mirror 
images over two mid-planes (two direct mirror images and the third is symmetric over the line 
of intersection of these two mid-planes), and these four points are associated to the same shape 
function. Consequently, the eight vertices of the cube are associated to the same shape function 
through symmetry. See Figure [2] for example. 



Our first tests are for Equation (|26j) with 5 = 4.0 and L = 0. According to Theorem 1.2 the 
optimal rate of convergence for the Finite Element solution should be obtained on triangulations 
with any k < 0.5, since rj — \fl/& + 4 > 1. The convergence rates e associated to triangulations 
with different values of k are listed in Table [T] Starting from an initial triangulation, we compute 
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Figure 2. The initial mesh on the unit cube (left); the mesh after one k refine- 
ment for the origin, k = 0.2 (right). 



j\e 


k = 0.1 


k = 0.2 


k = 0.3 


k = 0.4 


k = 0.5 


2 


0.42 


0.44 


0.56 


0.33 


-0.20 


3 


0.48 


0.68 


0.75 


0.79 


0.70 


4 


0.78 


0.81 


0.86 


0.88 


0.85 


5 


0.91 


0.92 


0.94 


0.95 


0.93 


6 


0.97 


0.97 


0.98 


0.99 


0.98 



Table 1 . Convergence rates e of finite element solutions solving equation ( 26 1 
with 8 = 4.0 and L = on different graded tetrahedra. 



the rates based on the comparison of the numerical errors on triangulations with consecutive 
^-refinements, 

row i \ v i-i ~ V M 
(27) e := log 2 , — l - , 

where Vj is the finite element solution on the mesh after j ^^refinements. Recall the dimension 
of the finite element space grows by a factor of 8 with one /c-refinement. By Theorem |1.2| for a 
sequence of optimal meshes, the error \v — Vj\^i is reduced by a factor of 2 for linear finite element 
approximations with each fc-refinement.Thus, e — > 1 implies that the optimal rate of convergence 
in Theorem 11.21 is achieved. 

Table [l] clearly shows that the convergence rates e approach 1 for all values of the grading 
parameter k. This is in agreement with our theory that the optimal rates of convergence are 
obtained for any triangulations with k < 0.5, since the singularity in the solution is not strong 
enough to be detectable for linear finite elements. 



In the second test, we implemented our method solving equation (26 1 with 5 = 0.6, L = 
and summarize the results in Table [i] Based on the upper bound r) = + 0.6 given in 

Theorem 1 1.2 1 we expect the optimal rate of convergence for the numerical solution as long as the 
grading parameter k < 2~ 1 /' ? S3 0.47. The convergence rates in Table [2] tend to 1 when k < 0.4, 
which implies the optimality of our finite element approximation on these meshes. However, 
when k — 0.5, the convergence rate is far less than 1 and there is a large gap between the rates 
corresponding to k = 0.4 and k = 0.5. This further confirms our theory that the upper bound of 
the suitable range of k for an optimal finite element approximation lies in (0.4,0.5). 



The third tests are for negative potentials in equation (26), where we set 8 = —0.1 and L = 20 



to satisfy the positivity requirement in Theorem 2.1 Our theoretical results indicate that the 



singularity in the solution due to the singular potential is stronger in this case and the optimal 

rate can be achieved only if the grading parameter k < 2 ss 0.167. Because of 

the limitation of the computation power, we only display the convergence results up to the 7th 
refinement for various graded parameters k in Table [3j We, however, still see the trend that 
appropriate gradings improve the convergence rate as predicted in Theorem |1.2| When k is close 
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k = 0.1 


k = 0.2 


k = 0.3 


k = 0.4 


k = 0.5 


2 


0.20 


0.30 


0.33 


0.11 


-0.03 


3 


0.54 


0.66 


0.69 


0.61 


0.39 


4 


0.74 


0.81 


0.83 


0.77 


0.60 


5 


0.88 


0.91 


0.92 


0.87 


0.72 


6 


0.95 


0.97 


0.98 


0.92 


0.79 



Table 2. Convergence rates e of finite element solutions solving equation (26) 
with S = 0.6 and L = on different graded tetrahedra. 



i\e 


k = 0.1 


k = 0.2 


k = 0.3 


k = 0.4 


k = 0.5 


2 


-0.10 


-0.05 


-0.09 


-0.16 


-0.03 


3 


0.32 


0.37 


0.30 


0.19 


0.07 


4 


0.51 


0.52 


0.44 


0.32 


0.18 


5 


0.67 


0.64 


0.53 


0.40 


0.26 


6 


0.80 


0.72 


0.59 


0.45 


0.32 



Table 3. Convergence rates e of finite element solutions solving equation (26) 
with S = —0.1 and L = 20 on different graded tetrahedra. 



to the optimal value 0.167 (i.e., k = 0.1 and 0.2), we have remarkable improvements. In particular, 
for k — 0.1, based on Table |3j we expect that the optimal rate occurs with further refinements. 
We have also implemented the method on graded meshes for the eigenvalue problem associated 



with equation (26), especially on the computation of the first eigenvalues. Namely, 

F k u := (-A + 5*(jr~ 2 )u = Xm 

on the unit cube, where Ai is the first eigenvalue of the operator. Depending on the choice of 6, 
the convergence rates for the numerical eigenvalues on graded meshes are roughly twice the rates 
for the numerical solutions of equation ( 26 ) (see Tables [l] [2j and |3| , and present similar trends 



for different gradings. 

All our numerical tests (Tables 1|2|3 and corresponding eigenvalue computations) convincingly 



verify Theorem |1.1| by comparing the rates of convergence for different singular potentials on 



different graded triangulations for the model operator in ( 26 ). The theoretical upper bounds 2 -1 '*' 
of the optimal range for the grading parameter k are also clearly demonstrated in these numerical 
results. In these tests, the initial triangulation of the unit cube consists of 12 tetrahedra and we 
consecutively refine the mesh using the fc-refinements up to level 7 that includes 12 x 8 7 ~ 2.5 x 10 7 
tetrahedra and roughly 4.2 million unknowns. Numerical experiments show that the condition 
numbers of our discrete systems grow by a factor of 4 for consecutive refinements, regardless of 
the value of k, which resembles the estimates given in [7] for the Laplace operator. However, the 
values of k affect the magnitude of the condition numbers. In general, smaller k leads to bad 
shapes for the tetrahedra and therefore results in larger condition numbers. The preconditioned 
conjugate gradient (PCG) method was used as the numerical solver for the discrete systems. 
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